Data Wrangling and Statistical Analysis of Assessment of Abundance, Aggregation, and Web Placement of Agelenopsis pennsylvanica
Brandi Pessman
2022-11-20
Setting the Working Directory and Global Code Chunk Options
Using include = FALSE means that this chunk will run, but the code and the results will not appear in the knitted document. We have set message and warning to FALSE for every chunk so the knitted document is not bombarded with messages and warnings.
Packages to Load
There is a brief explanation of each package’s purpose left as a comment following the library call.
Import Data
The document “webs” is a text file in the data folder of the repository/directory.
Data Wrangling
Here, we rename some levels of factors Location and Land and put them in the order we want the categories to appear in a graph. We also create another data set “sites” that restricts “webs” to just a single data point for each site. We use “sites” to compare search distance, the number of webs, and the number of spiders. We will use “webs” for web height and nearest neighbor analyses.
At one time, we tried scaling and centering the data in difference ways. I want to make mention of it here because it might be useful information for me to look back on in the future. Scaling is dividing each datum by the standard deviation of the data and centering is subtracting the mean of the data from each datum to try to normalize the data; a second type of centering and scaling involves subtracting the median and dividing by the interquartile range; need to remember to back transform when making figures; we didn’t scale and center because log-transformations on zero-inflated variables were enough to remove convergence errors
Transformation of Variables
Many data are zero-inflated or right skewed. We check the distribution of the data and check distributions after log-transformation for right-skewed data not containing zeros. Echo = FALSE shows the results, but does not show the code used.
##
## Shapiro-Wilk normality test
##
## data: sites$tree100m
## W = 0.79926, p-value = 0.0004871
##
## Shapiro-Wilk normality test
##
## data: sites$imperv100m
## W = 0.74128, p-value = 6.883e-05
##
## Shapiro-Wilk normality test
##
## data: sites$TotalSub
## W = 0.85215, p-value = 0.003707##
## Shapiro-Wilk normality test
##
## data: log(sites$TotalSub)
## W = 0.96252, p-value = 0.5417
##
## Shapiro-Wilk normality test
##
## data: sites$Traffic_Dist_Road
## W = 0.70845, p-value = 2.511e-05##
## Shapiro-Wilk normality test
##
## data: log(sites$Traffic_Dist_Road)
## W = 0.93619, p-value = 0.1651
##
## Shapiro-Wilk normality test
##
## data: sites$Traffic_Dist_Highway
## W = 0.59575, p-value = 1.195e-06##
## Shapiro-Wilk normality test
##
## data: log(sites$Traffic_Dist_Highway)
## W = 0.9091, p-value = 0.04539
##
## Shapiro-Wilk normality test
##
## data: sites$spec_rad
## W = 0.85174, p-value = 0.003645##
## Shapiro-Wilk normality test
##
## data: log(sites$spec_rad)
## W = 0.84964, p-value = 0.003345
##
## Shapiro-Wilk normality test
##
## data: sites$light_rad
## W = 0.84052, p-value = 0.002319##
## Shapiro-Wilk normality test
##
## data: log(sites$light_rad)
## W = 0.81646, p-value = 0.0009161
##
## Shapiro-Wilk normality test
##
## data: sites$patch_area_mm
## W = 0.70203, p-value = 2.077e-05##
## Shapiro-Wilk normality test
##
## data: log(sites$patch_area_mm)
## W = 0.89878, p-value = 0.0281
##
## Shapiro-Wilk normality test
##
## data: sites$road_length_m
## W = 0.87899, p-value = 0.01156
##
## Shapiro-Wilk normality test
##
## data: sites$trail_length_m
## W = 0.7756, p-value = 0.0002127We conclude to use log-transformations for TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, and patch_area_km.
Our analyses include looking at variation within each location. As such, this chunk subsets the “sites” and “webs” data by Location.
Variable Reduction Based on Correlations
Since many of these variables are likely highly correlated, we want to reduce the variables that we include in our model by removing variables that are highly correlated.
corr <- sites %>%
select(tree100m, imperv100m, TotalSub,
Traffic_Dist_Road, Traffic_Dist_Highway,
spec_rad, light_rad, patch_area_mm, road_length_m)
# We exclude trail length from global correlations because trail length can only collected in Wilderness Park, and not for UNL City Campus
# When running findCorrelations, we get a vector of variables to remove to reduce pairwise correlations.
findCorrelation(cor(corr, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 6 and column 2 with corr 0.73
## Means: 0.693 vs 0.594 so flagging column 6
## Compare row 2 and column 8 with corr 0.77
## Means: 0.68 vs 0.568 so flagging column 2
## Compare row 8 and column 9 with corr 0.738
## Means: 0.633 vs 0.53 so flagging column 8
## Compare row 9 and column 7 with corr 0.706
## Means: 0.613 vs 0.487 so flagging column 9
## Compare row 7 and column 1 with corr 0.623
## Means: 0.601 vs 0.44 so flagging column 7
## Compare row 1 and column 3 with corr 0.662
## Means: 0.494 vs 0.336 so flagging column 1
## All correlations <= 0.6## [1] "spec_rad" "imperv100m" "patch_area_mm" "road_length_m"
## [5] "light_rad" "tree100m"# This suggests that we remove tree100m, imperv100, spec_rad, light_rad, patch_area_km, road_length_m
# That leaves us with TotalSub, Traffic_Dist_Road, and Traffic_Dist_Highway
chart.Correlation(corr, histogram = TRUE, method = "spearman") # kendall, spearman
corr <- sites %>%
select(TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway)
chart.Correlation(corr, histogram = TRUE, method = "spearman") # kendall, spearman
# Let's also test the variables after transformation
corr_transformed <- sites %>%
select(tree100m, imperv100m, log_TotalSub,
log_tdr, log_tdh,
spec_rad, light_rad, log_patch, road_length_m)
findCorrelation(cor(corr_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 6 and column 2 with corr 0.73
## Means: 0.693 vs 0.594 so flagging column 6
## Compare row 2 and column 8 with corr 0.77
## Means: 0.68 vs 0.568 so flagging column 2
## Compare row 8 and column 9 with corr 0.738
## Means: 0.633 vs 0.53 so flagging column 8
## Compare row 9 and column 7 with corr 0.706
## Means: 0.613 vs 0.487 so flagging column 9
## Compare row 7 and column 1 with corr 0.623
## Means: 0.601 vs 0.44 so flagging column 7
## Compare row 1 and column 3 with corr 0.662
## Means: 0.494 vs 0.336 so flagging column 1
## All correlations <= 0.6## [1] "spec_rad" "imperv100m" "log_patch" "road_length_m"
## [5] "light_rad" "tree100m"# This suggests to remove tree100m, imperv100m, spec_rad, light_rad, log_patch, and road_length_m.
# This leaves us with log_TotalSub, log_tdr, and log_tdh, the same results as the non-transformed.
chart.Correlation(corr_transformed, histogram = TRUE, method = "spearman") # kendall, spearman
corr_transformed <- sites %>%
select(log_TotalSub, log_tdr, log_tdh)
chart.Correlation(corr_transformed, histogram = TRUE, method = "spearman") # kendall, spearman
corr_forest <- sites %>%
filter(Location == "Urban Forest") %>%
select(tree100m, imperv100m, TotalSub,
Traffic_Dist_Road, Traffic_Dist_Highway,
spec_rad, light_rad, patch_area_km, road_length_m, trail_length_m)
# Notice trail length is included
findCorrelation(cor(corr_forest, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 10 and column 2 with corr 0.696
## Means: 0.443 vs 0.315 so flagging column 10
## Compare row 2 and column 9 with corr 0.899
## Means: 0.361 vs 0.286 so flagging column 2
## Compare row 9 and column 1 with corr 0.679
## Means: 0.279 vs 0.272 so flagging column 9
## Compare row 6 and column 7 with corr 0.685
## Means: 0.371 vs 0.258 so flagging column 6
## Compare row 7 and column 3 with corr 0.911
## Means: 0.296 vs 0.219 so flagging column 7
## All correlations <= 0.6## [1] "trail_length_m" "imperv100m" "road_length_m" "spec_rad"
## [5] "light_rad"# This suggests removing imperv100m, spec_rad, light_rad, road_length_m, and trail_length_m
# We will keep tree100m, TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, and patch_area_km
chart.Correlation(corr_forest, histogram = TRUE, method = "spearman")
corr_forest <- sites %>%
filter(Location == "Urban Forest") %>%
select(tree100m, TotalSub, Traffic_Dist_Road, patch_area_km, Traffic_Dist_Highway)
chart.Correlation(corr_forest, histogram = TRUE, method = "spearman")
# Let's try with the log-transformed
corr_forest_transformed <- sites %>%
filter(Location == "Urban Forest") %>%
select(tree100m, imperv100m, log_TotalSub,
log_tdr, log_tdh,
spec_rad, light_rad, log_patch, road_length_m, trail_length_m)
findCorrelation(cor(corr_forest_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 10 and column 2 with corr 0.696
## Means: 0.443 vs 0.315 so flagging column 10
## Compare row 2 and column 9 with corr 0.899
## Means: 0.361 vs 0.286 so flagging column 2
## Compare row 9 and column 1 with corr 0.679
## Means: 0.279 vs 0.272 so flagging column 9
## Compare row 6 and column 7 with corr 0.685
## Means: 0.371 vs 0.258 so flagging column 6
## Compare row 7 and column 3 with corr 0.911
## Means: 0.296 vs 0.219 so flagging column 7
## All correlations <= 0.6## [1] "trail_length_m" "imperv100m" "road_length_m" "spec_rad"
## [5] "light_rad"# The same variables are dropped and kept
corr_forest_transformed <- sites %>%
filter(Location == "Urban Forest") %>%
select(tree100m, log_TotalSub, log_tdr, patch_area_km, log_tdh)
chart.Correlation(corr_forest_transformed, histogram = TRUE, method = "spearman")
corr_center <- sites %>%
filter(Location == "Urban Center") %>%
select(tree100m, imperv100m, TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, spec_rad, light_rad, patch_area_km, road_length_m)
# Trail length is removed because we could not measure trail length on campus
findCorrelation(cor(corr_center, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 1 and column 2 with corr 0.754
## Means: 0.477 vs 0.322 so flagging column 1
## Compare row 2 and column 4 with corr 0.636
## Means: 0.388 vs 0.282 so flagging column 2
## Compare row 7 and column 5 with corr 0.951
## Means: 0.352 vs 0.252 so flagging column 7
## All correlations <= 0.6## [1] "tree100m" "imperv100m" "light_rad"# This suggests that we remove tree100m, imperv100m, and light_rad
# We will keep TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, spec_rad, patch_area_km, and road_length_m
chart.Correlation(corr_center, histogram = TRUE, method = "spearman")
# Let's try with the transformed
corr_center_transformed <- sites %>%
filter(Location == "Urban Center") %>%
select(tree100m, imperv100m, log_TotalSub, log_tdr, log_tdh, spec_rad, light_rad, log_patch, road_length_m)
findCorrelation(cor(corr_center_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) ## Compare row 1 and column 2 with corr 0.754
## Means: 0.477 vs 0.322 so flagging column 1
## Compare row 2 and column 4 with corr 0.636
## Means: 0.388 vs 0.282 so flagging column 2
## Compare row 7 and column 5 with corr 0.951
## Means: 0.352 vs 0.252 so flagging column 7
## All correlations <= 0.6## [1] "tree100m" "imperv100m" "light_rad"chart.Correlation(corr_center_transformed, histogram = TRUE, method = "pearson")
# The same variables are kept or removedFor overall analysis of predictors, we will include:
log_TotalSub
log_tdr
log_tdh
For the urban forest subset, we will use:
tree100m
log_TotalSub
log_tdr
log_tdh
log_patch
Finally, for the urban center subset, we will use:
log_TotalSub
log_tdr
log_tdh
spec_rad
log_patch
road_length_m
Graphs, Stats, and Assumptions Comparing Predictors by Location
Here, we just want to see if each predictor differs between the two Locations: urban forest and urban center. We build models, check assumptions of the model, and make graphs to represent differences.
##
## Call:
## glm(formula = tree100m ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4483 -1.7389 -0.8778 0.9780 2.8970
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.85340 0.04605 83.68 <2e-16 ***
## LocationUrban Center -3.39388 0.23399 -14.50 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 664.364 on 21 degrees of freedom
## Residual deviance: 58.556 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 6## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 47.2 3.26
## 2 Urban Center 1.58 0.692
##
## Call:
## glm(formula = imperv100m ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.92734 -1.59531 0.02675 0.71981 2.67358
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2410 0.2803 0.86 0.39
## LocationUrban Center 4.0965 0.2823 14.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1039.415 on 21 degrees of freedom
## Residual deviance: 40.746 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 6## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 1.27 0.644
## 2 Urban Center 76.5 2.70
##
## Call:
## glm(formula = log_tdr ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.07303 -0.43708 -0.08942 0.38719 1.04170
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3786 0.1587 8.685 <2e-16 ***
## LocationUrban Center 0.2841 0.2025 1.403 0.161
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8.5848 on 21 degrees of freedom
## Residual deviance: 6.5798 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 3.97 0.349
## 2 Urban Center 5.27 0.393
##
## Call:
## glm(formula = log_tdh ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.41615 -0.19952 -0.00187 0.12375 0.69255
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2585 0.1685 7.467 8.22e-14 ***
## LocationUrban Center 0.1178 0.2224 0.530 0.596
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1.6528 on 21 degrees of freedom
## Residual deviance: 1.3707 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 3.52 0.116
## 2 Urban Center 3.96 0.179
##
## Call:
## glm(formula = log_TotalSub ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.47861 -0.36634 0.02223 0.26908 0.61117
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9173 0.1999 4.589 4.45e-06 ***
## LocationUrban Center -0.8283 0.3409 -2.430 0.0151 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 10.6306 on 21 degrees of freedom
## Residual deviance: 4.3622 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 2.50 0.139
## 2 Urban Center 1.09 0.148
##
## Call:
## glm(formula = spec_rad ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.73540 -0.10657 -0.01758 0.09970 0.75755
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.92321 0.02697 182.520 < 2e-16 ***
## LocationUrban Center 0.10975 0.03565 3.079 0.00208 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 11.0025 on 21 degrees of freedom
## Residual deviance: 1.4789 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 3## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 137. 0.434
## 2 Urban Center 153. 1.26
##
## Call:
## glm(formula = light_rad ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.1501 -1.0590 -0.5849 1.8735 4.9674
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7725 0.1303 13.60 <2e-16 ***
## LocationUrban Center 2.9888 0.1331 22.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1421.64 on 21 degrees of freedom
## Residual deviance: 121.73 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 5.89 0.956
## 2 Urban Center 117. 9.80
##
## Call:
## glm(formula = log_patch ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.74174 -0.19285 0.01795 0.11564 0.56118
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.35374 0.09747 24.147 < 2e-16 ***
## LocationUrban Center -0.54815 0.15231 -3.599 0.00032 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 15.4711 on 21 degrees of freedom
## Residual deviance: 2.2471 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 10.5 0.224
## 2 Urban Center 6.08 0.286
##
## Call:
## glm(formula = road_length_m ~ Location, family = "poisson", data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.585 -5.585 -1.065 2.315 8.101
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.74702 0.08007 34.31 <2e-16 ***
## LocationUrban Center 1.78585 0.08548 20.89 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1098.64 on 21 degrees of freedom
## Residual deviance: 446.12 on 20 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 6## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 15.6 7.97
## 2 Urban Center 93.0 6.93
Response Variable: Search Distance
Search Distance by Location
ggarrange(ggplot(sites, aes(x = WalkDist)) +
geom_density(),
ggplot(sites, aes(x = log(WalkDist))) +
geom_density(), ncol = 2)
shapiro.test(sites$WalkDist) # not normal##
## Shapiro-Wilk normality test
##
## data: sites$WalkDist
## W = 0.86946, p-value = 0.007649shapiro.test(log(sites$WalkDist)) # normal##
## Shapiro-Wilk normality test
##
## data: log(sites$WalkDist)
## W = 0.94443, p-value = 0.244search <- glm(WalkDist ~ Location, data = sites, family = poisson)
summary(search)##
## Call:
## glm(formula = WalkDist ~ Location, family = poisson, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.092 -7.441 -2.375 3.649 12.920
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.69410 0.03025 155.2 <2e-16 ***
## LocationUrban Center -0.46605 0.04615 -10.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1489.7 on 21 degrees of freedom
## Residual deviance: 1386.4 on 20 degrees of freedom
## AIC: 1517.2
##
## Number of Fisher Scoring iterations: 5search2 <- glm(log(WalkDist) ~ Location, data = sites)
summary(search2)##
## Call:
## glm(formula = log(WalkDist) ~ Location, data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2849 -0.6105 0.2763 0.9588 1.6204
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.9780 0.4062 9.792 4.5e-09 ***
## LocationUrban Center -0.1102 0.5501 -0.200 0.843
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.650385)
##
## Null deviance: 33.074 on 21 degrees of freedom
## Residual deviance: 33.008 on 20 degrees of freedom
## AIC: 77.359
##
## Number of Fisher Scoring iterations: 2sites %>%
group_by(Location) %>%
summarize(mean = mean(WalkDist),
se = plotrix::std.error(WalkDist))## # A tibble: 2 × 3
## Location mean se
## <fct> <dbl> <dbl>
## 1 Urban Forest 109. 31.7
## 2 Urban Center 68.6 16.2
Search Distance by Predictors
# This is one way to look at the overall data, but I think I end up not using this information
models = list(WalkDist ~ 1,
WalkDist ~ log_TotalSub,
WalkDist ~ log_tdr,
WalkDist ~ log_tdh,
WalkDist ~ log_TotalSub + log_tdr,
WalkDist ~ log_TotalSub + log_tdh,
WalkDist ~ log_tdr + log_tdh,
WalkDist ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = sites, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web Overall.", split.tables = Inf)| Modnames | K | AICc | Delta_AICc | ModelLik | AICcWt | LL | Cum.Wt | |
|---|---|---|---|---|---|---|---|---|
| 8 | log_TotalSub + log_tdr + log_tdh | 4 | 1504 | 0 | 1 | 0.9999 | -747.1 | 0.9999 |
| 5 | log_TotalSub + log_tdr | 3 | 1523 | 18.94 | 7.709e-05 | 7.709e-05 | -758.1 | 1 |
| 6 | log_TotalSub + log_tdh | 3 | 1558 | 53.06 | 3.013e-12 | 3.013e-12 | -775.1 | 1 |
| 2 | log_TotalSub | 2 | 1565 | 60.03 | 9.222e-14 | 9.221e-14 | -779.9 | 1 |
| 3 | log_tdr | 2 | 1597 | 92.5 | 8.208e-21 | 8.207e-21 | -796.2 | 1 |
| 7 | log_tdr + log_tdh | 3 | 1600 | 95.14 | 2.187e-21 | 2.187e-21 | -796.2 | 1 |
| 1 | 1 | 1 | 1619 | 114.2 | 1.616e-25 | 1.616e-25 | -808.2 | 1 |
| 4 | log_tdh | 2 | 1621 | 116.6 | 4.912e-26 | 4.911e-26 | -808.2 | 1 |
# The top model keeps all 3 variables and no other model is within 2 Model Likelihood Units
summary(glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson))##
## Call:
## glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh, family = poisson,
## data = sites)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -14.463 -7.694 -1.934 5.088 15.497
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.45699 0.24451 10.049 < 2e-16 ***
## log_TotalSub 0.31671 0.03252 9.740 < 2e-16 ***
## log_tdr 0.13279 0.01775 7.480 7.42e-14 ***
## log_tdh 0.21434 0.04469 4.796 1.62e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1489.7 on 21 degrees of freedom
## Residual deviance: 1367.3 on 18 degrees of freedom
## AIC: 1502.1
##
## Number of Fisher Scoring iterations: 5global <- glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = F, evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_o)
info <- out[[12]]
#TotalSub (log_TotalSub)
exp(info[4, 1])## [1] 1.372605exp(info[4, 2])## [1] 1.033052#Traffic_Dist_Road (log_tdr)
exp(info[3, 1])## [1] 1.142011exp(info[3, 2])## [1] 1.017911#Traffic_Dist_Highway (log_tdh)
exp(info[2, 1])## [1] 1.239044exp(info[2, 2])## [1] 1.045706with(summary(top_model_wd_o), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.0821148#summary(model.avg(model_dredge, subset = delta <= 2))
#summary(model.avg(model_dredge))
model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh, family = poisson,
## data = sites)
## ---
## Model selection table
## (Int) log_tdh log_tdr log_TtS df logLik AICc delta weight
## 8 2.457 0.214300 0.13280 0.3167 4 -747.072 1504.5 0.00 1
## 7 3.489 0.11210 0.2462 3 -758.052 1523.4 18.94 0
## 6 3.488 0.142400 0.2455 3 -775.110 1557.6 53.06 0
## 5 4.094 0.2063 2 -779.947 1564.5 60.03 0
## 3 4.079 0.08155 2 -796.182 1597.0 92.50 0
## 4 4.040 0.009906 0.08199 3 -796.153 1599.6 95.14 0
## 1 4.467 1 -808.233 1618.7 114.17 0
## 2 4.502 -0.009444 2 -808.208 1621.0 116.55 0
## Models ranked by AICc(x)car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2## log_tdh log_tdr log_TotalSub
## 1.301594 1.124043 1.357283# urban forest
global <- glm(WalkDist ~ tree100m + log_TotalSub + log_tdh + log_patch + log_tdr, data = sites_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = F, evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_f)
info <- out[[12]]
#tree100m
info[5, 1]## [1] 0.04991404info[5, 2]## [1] 0.004240799#TotalSub
-exp(info[4, 1])## [1] -0.3642043exp(info[4, 2])## [1] 1.099856#Traffic_Dist_Road
exp(info[3, 1])## [1] 2.051607exp(info[3, 2])## [1] 1.039132#Traffic_Dist_Highway
exp(info[2, 1])## [1] 3.875139exp(info[2, 2])## [1] 1.116225with(summary(top_model_wd_f), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.5328974model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = WalkDist ~ tree100m + log_TotalSub + log_tdh +
## log_patch + log_tdr, family = poisson, data = sites_forest)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS t10 df logLik AICc delta
## 31 -2.9620 1.3550 0.7186 -1.0100 0.049910 5 -237.637 500.3 0.00
## 32 -4.2870 0.111800 1.3660 0.7309 -0.9987 0.050600 6 -235.218 510.4 10.16
## 23 -4.9080 1.4220 0.5793 0.045490 4 -298.101 612.2 111.93
## 24 -6.3640 0.121500 1.4390 0.5930 0.046770 5 -295.358 615.7 115.44
## 15 1.3710 0.9436 0.4754 -0.8079 4 -311.229 638.5 138.18
## 16 -0.1680 0.130700 0.9870 0.4797 -0.8106 5 -308.126 641.3 140.98
## 29 3.7220 0.6233 -1.1830 0.027800 4 -332.487 681.0 180.70
## 30 3.4360 0.025250 0.6257 -1.1830 0.028000 5 -332.338 689.7 189.40
## 7 -0.7565 1.1140 0.3608 3 -355.948 721.9 221.62
## 8 -2.0680 0.111200 1.1500 0.3635 4 -353.705 723.4 223.14
## 13 5.3650 0.4481 -1.0220 3 -368.601 747.2 246.93
## 14 5.1990 0.015970 0.4485 -1.0240 4 -368.545 753.1 252.82
## 21 2.0920 0.3613 0.023610 3 -413.559 837.1 336.84
## 22 1.9760 0.010040 0.3621 0.023740 4 -413.537 843.1 342.80
## 27 1.6500 0.9880 -0.3968 0.010370 4 -414.777 845.6 345.28
## 11 2.4340 0.8922 -0.3759 3 -418.557 847.1 346.84
## 12 2.0780 0.029380 0.9042 -0.3742 4 -418.359 852.7 352.44
## 28 1.7170 -0.006520 0.9872 -0.3978 0.010520 5 -414.767 854.5 354.26
## 19 0.7985 0.9742 0.008801 3 -428.397 866.8 366.52
## 3 1.4210 0.9170 2 -431.230 868.2 367.90
## 4 1.0170 0.034130 0.9295 3 -430.969 871.9 371.66
## 20 0.7048 0.009042 0.9759 0.008645 4 -428.379 872.8 372.48
## 5 3.7340 0.2338 2 -437.478 880.7 380.40
## 6 3.9550 -0.021130 0.2341 3 -437.380 884.8 384.49
## 9 5.7640 -0.4343 2 -457.538 920.8 420.52
## 10 6.1830 -0.039100 -0.4371 3 -457.142 924.3 424.01
## 25 5.6580 -0.4322 0.002144 3 -457.300 924.6 424.33
## 26 6.1230 -0.046060 -0.4359 0.002736 4 -456.763 929.5 429.25
## 1 4.6940 1 -475.501 953.5 453.23
## 17 4.5250 0.003571 2 -474.832 955.4 455.10
## 2 5.0660 -0.035390 2 -475.193 956.1 455.83
## 18 4.9430 -0.041030 0.003855 3 -474.415 958.8 458.56
## weight
## 31 0.994
## 32 0.006
## 23 0.000
## 24 0.000
## 15 0.000
## 16 0.000
## 29 0.000
## 30 0.000
## 7 0.000
## 8 0.000
## 13 0.000
## 14 0.000
## 21 0.000
## 22 0.000
## 27 0.000
## 11 0.000
## 12 0.000
## 28 0.000
## 19 0.000
## 3 0.000
## 4 0.000
## 20 0.000
## 5 0.000
## 6 0.000
## 9 0.000
## 10 0.000
## 25 0.000
## 26 0.000
## 1 0.000
## 17 0.000
## 2 0.000
## 18 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2## log_tdh log_tdr log_TotalSub tree100m
## 1.444130 1.612738 1.259335 1.840585# urban center
global <- glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = sites_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = F, evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_c)
info <- out[[12]]
#TotalSub
exp(info[5, 1])## [1] 1.726exp(info[5, 2])## [1] 1.084755#patch_area_km
exp(info[2, 1])## [1] 2.864928exp(info[2, 2])## [1] 1.060802#spec_rad
info[6, 1]## [1] 0.1246571info[6, 2]## [1] 0.01458074#Traffic_Dist_Road
exp(info[4, 1])## [1] 1.573514exp(info[4, 2])## [1] 1.044837#Traffic_Dist_Highway
exp(info[3, 1])## [1] 1.786049exp(info[3, 2])## [1] 1.095576with(summary(top_model_wd_c), 1 - deviance/null.deviance) # This gives the R squared value## [1] 0.9399876model.sel(model_dredge) #estimates same signs## Global model call: glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh + spec_rad +
## log_patch + road_length_m, family = poisson, data = sites_center)
## ---
## Model selection table
## (Int) log_ptc log_tdh log_tdr log_TtS rod_lng_m spc_rad df logLik
## 48 -26.91000 1.0530 0.58000 0.4533 0.545800 0.124700 6 -49.077
## 64 -27.07000 1.0540 0.58820 0.4582 0.544400 4.438e-04 0.125100 7 -49.045
## 46 -18.79000 0.9700 0.3030 0.467300 0.096030 5 -69.398
## 62 -18.05000 0.9688 0.2870 0.477700 -2.947e-03 0.093480 6 -67.840
## 40 -23.31000 0.9764 0.50860 0.5486 0.106700 5 -72.569
## 56 -23.55000 0.9785 0.53010 0.5595 9.858e-04 0.106600 6 -72.413
## 38 -16.35000 0.9141 0.4031 0.082200 4 -87.897
## 54 -16.16000 0.9160 0.3926 -2.429e-03 0.082660 5 -86.790
## 16 -4.79700 0.8060 0.35770 0.3880 0.347000 5 -92.534
## 32 -4.18200 0.7949 0.30640 0.3545 0.376000 -2.056e-03 6 -91.779
## 30 -2.16500 0.7837 0.2588 0.371000 -4.289e-03 5 -97.463
## 14 -2.78000 0.7983 0.2954 0.306900 4 -101.490
## 8 -4.47000 0.7813 0.32270 0.4521 4 -102.558
## 24 -4.43700 0.7806 0.31990 0.4506 -9.962e-05 5 -102.556
## 6 -2.69600 0.7794 0.3634 3 -109.821
## 58 -11.68000 0.8625 0.840500 -6.534e-03 0.066040 5 -104.636
## 22 -2.30800 0.7684 0.3491 -2.607e-03 4 -108.248
## 60 -11.77000 0.8623 0.01075 0.845300 -6.509e-03 0.066310 6 -104.625
## 42 -13.35000 0.8539 0.825900 0.073500 4 -114.026
## 44 -13.62000 0.8523 0.03948 0.841100 0.074180 5 -113.881
## 26 -0.33890 0.7077 0.716900 -7.737e-03 4 -124.349
## 28 -0.04446 0.7153 -0.07423 0.688600 -7.937e-03 5 -123.800
## 10 -0.84410 0.6920 0.625600 3 -140.065
## 12 -0.73870 0.6951 -0.02912 0.616400 4 -139.983
## 52 -2.37500 0.6594 -0.24640 -6.566e-03 0.026110 5 -163.515
## 20 1.97600 0.6132 -0.25760 -6.713e-03 4 -166.678
## 50 -3.52000 0.6142 -5.326e-03 0.028380 4 -169.190
## 18 1.17400 0.5606 -5.367e-03 3 -173.031
## 36 -3.50600 0.6315 -0.17530 0.028930 4 -173.540
## 34 -4.11800 0.6026 0.029610 3 -176.457
## 4 1.21900 0.5873 -0.17900 3 -177.822
## 2 0.72240 0.5547 2 -180.946
## 21 3.78000 0.1896 -6.382e-03 3 -246.767
## 53 6.51000 0.1793 -6.883e-03 -0.017160 4 -245.220
## 29 3.80300 0.1709 0.122100 -7.024e-03 4 -245.751
## 23 3.34400 0.07582 0.2052 -5.828e-03 4 -246.089
## 61 6.43400 0.1612 0.117700 -7.549e-03 -0.016480 5 -244.290
## 55 5.95200 0.06141 0.1931 -6.387e-03 -0.015890 5 -244.790
## 31 3.42200 0.06619 0.1858 0.113100 -6.496e-03 5 -245.228
## 7 2.43700 0.14180 0.2262 3 -252.264
## 5 3.15300 0.1976 2 -254.778
## 63 5.96200 0.05205 0.1739 0.110600 -7.086e-03 -0.015420 6 -243.978
## 39 3.39400 0.14100 0.2231 -0.006122 4 -252.053
## 37 4.17700 0.1940 -0.006552 3 -254.531
## 15 2.43500 0.14110 0.2235 0.017190 4 -252.240
## 13 3.14700 0.1940 0.023060 3 -254.735
## 45 4.12500 0.1922 0.013050 -0.006248 4 -254.518
## 47 3.36500 0.14070 0.2218 0.008304 -0.005939 5 -252.048
## 57 8.67200 0.317300 -9.203e-03 -0.025840 4 -257.913
## 59 9.27400 -0.09140 0.304500 -9.818e-03 -0.026950 5 -256.705
## 25 4.59900 0.343500 -8.289e-03 3 -262.339
## 27 4.96100 -0.07796 0.334600 -8.768e-03 4 -261.447
## 49 9.56200 -7.437e-03 -0.030380 3 -265.893
## 51 10.32000 -0.11640 -8.312e-03 -0.031830 4 -264.034
## 17 4.82400 -6.523e-03 2 -271.928
## 19 5.27900 -0.09958 -7.195e-03 3 -270.558
## 9 3.96600 0.233400 2 -275.865
## 41 6.01800 0.202000 -0.013160 3 -274.663
## 11 4.02500 -0.01410 0.231200 3 -275.836
## 33 7.22500 -0.019560 2 -278.356
## 43 6.06900 -0.01300 0.200000 -0.013140 4 -274.639
## 1 4.22800 1 -281.081
## 35 7.33500 -0.03213 -0.019450 3 -278.209
## 3 4.36800 -0.03536 2 -280.903
## AICc delta weight
## 48 127.0 0.00 0.999
## 64 140.1 13.14 0.001
## 46 158.8 31.84 0.000
## 62 164.5 37.53 0.000
## 40 165.1 38.18 0.000
## 56 173.6 46.67 0.000
## 38 189.5 62.55 0.000
## 54 193.6 66.63 0.000
## 16 205.1 78.11 0.000
## 32 212.4 85.40 0.000
## 30 214.9 87.97 0.000
## 14 216.7 89.74 0.000
## 8 218.8 91.88 0.000
## 24 225.1 98.16 0.000
## 6 228.6 101.69 0.000
## 58 229.3 102.32 0.000
## 22 230.2 103.26 0.000
## 60 238.0 111.10 0.000
## 42 241.8 114.81 0.000
## 44 247.8 120.81 0.000
## 26 262.4 135.46 0.000
## 28 267.6 140.65 0.000
## 10 289.1 162.18 0.000
## 12 293.7 166.73 0.000
## 52 347.0 220.08 0.000
## 20 347.1 220.12 0.000
## 50 352.1 225.14 0.000
## 18 355.1 228.11 0.000
## 36 360.8 233.84 0.000
## 34 361.9 234.96 0.000
## 4 364.6 237.69 0.000
## 2 367.2 240.27 0.000
## 21 502.5 375.58 0.000
## 53 504.2 377.20 0.000
## 29 505.2 378.26 0.000
## 23 505.9 378.94 0.000
## 61 508.6 381.63 0.000
## 55 509.6 382.63 0.000
## 31 510.5 383.50 0.000
## 7 513.5 386.57 0.000
## 5 514.9 387.94 0.000
## 63 516.8 389.80 0.000
## 39 517.8 390.87 0.000
## 37 518.1 391.11 0.000
## 15 518.2 391.24 0.000
## 13 518.5 391.52 0.000
## 45 522.8 395.80 0.000
## 47 524.1 397.14 0.000
## 57 529.5 402.59 0.000
## 59 533.4 406.46 0.000
## 25 533.7 406.72 0.000
## 27 536.6 409.65 0.000
## 49 540.8 413.83 0.000
## 51 541.8 414.83 0.000
## 17 549.2 422.24 0.000
## 19 550.1 423.16 0.000
## 9 557.1 430.11 0.000
## 41 558.3 431.37 0.000
## 11 560.7 433.72 0.000
## 33 562.0 435.09 0.000
## 43 563.0 436.04 0.000
## 1 564.6 437.61 0.000
## 35 565.4 438.46 0.000
## 3 567.1 440.19 0.000
## Models ranked by AICc(x)#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) # road and highway disturbance between 2 and 3## log_patch log_tdh log_tdr log_TotalSub spec_rad
## 1.874528 2.076776 2.677800 1.554569 1.661330